clear all
clc
muc=0.033;
C=[10:1:40];
Lambdac= 8;
T=10;
R=zeros(length(C),1);
numbers=zeros(length(C),1);

for j=1:length(C)
n=ceil(((Lambdac - muc)*T)/(T*C(j)*muc-1));
numbers(j)= n;
  Pv=decreasingdist(n);
%    Pv=Gaussiandist(n);
Pc=Gaussiandist(n);
lb=zeros(n,1);
ub=ones(n,1);
ub(n+1)=inf;
b1=zeros(n+2,1);
lb(n+1)=Lambdac + n*(1/T) ;
for i=1:n
    b1(i)= Lambdac*Pc(i) + (1/T) ;
end
b1(n+1)= -C(j)*n*muc;
b1(n+2)= -muc;
cond = @(x)parameterfun(x,n,Pv,b1);
obj=@(x)x(n+1);
x0 = ones(1,n);
x0(n+1)= Lambdac;
[x,fval,exitflag,output,lambda] = fmincon(obj,x0,[],[],[],[],lb,ub,cond);
R(j)= x(n+1);
[cout,ceqout] = parameterfun(x,n,Pv,b1);
% Waitings=zeros(n,1);
% for k=1:n
%    Waitings(k)=1/(-cout(k)+(1/T(j))) ;
% end
% WaitingsMax(j)=max(Waitings);
% WaitingsMin(j)=min(Waitings);
end


plot(C,R)
hold on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Using the max n for all of them 
% clear all
% clc
% muc=0.033;
% C=40;
% Lambdac= 5;
% T=2:0.1:15;
% n=ceil(((Lambdac - muc)*min(T))/(min(T)*C*muc-1));
% R=zeros(length(T),1);
% Waitings=zeros(length(T),1);
% % Pv=Gaussiandist(n);
% Pv=decreasingdist(n);
% Pc=Gaussiandist(n);
% lb=zeros(n,1);
% ub=ones(n,1);
% ub(n+1)=inf;
% b1=zeros(n+2,1);
% for j=1:length(T)
% lb(n+1)=Lambdac + n*(1/T(j)) ;
% for i=1:n
%     b1(i)= Lambdac*Pc(i) + (1/T(j)) ;
% end
% b1(n+1)= -C*n*muc;
% b1(n+2)= -muc;
% cond = @(x)parameterfun(x,n,Pv,b1);
% obj=@(x)x(n+1);
% x0 = ones(1,n);
% x0(n+1)= Lambdac;
% [x,fval,exitflag,output,lambda] = fmincon(obj,x0,[],[],[],[],lb,ub,cond);
% R(j)= x(n+1);
% [cout,ceqout] = parameterfun(x,n,Pv,b1);
% Waitings(j)=max(1/(-cout(1:n)+0.1));
% end
% 
% plot(T,R);
% hold on